
#include "ENSideElement.h"
#include "RayLine.h"
#include "LineSegment.h"
#include "MooseRandom.h"
#include "libmesh/fe_interface.h"
#include "libmesh/fe_type.h"

using namespace std;

ENSideElement::ENSideElement(const Elem *elem, const Point normal, int devide, Real transmissivity, Real absorptivity) :
	_elem(elem),
	_normal(normal),
	_devide(devide),
	_transmissivity(transmissivity),
	_absorptivity(absorptivity)
{
	Disintegrate();
}

void ENSideElement::Disintegrate()
{
	Point p = Point();
	unsigned int dim = _elem->dim();

	_Disintegrate_aera.clear();
	_Disintegrate_centre.clear();

	if (_devide < 1)
	{
		mooseError("Cannot disintegrate! ");
	}

	if (dim == 2) //3D
	{
		if (_elem->type() == TRI3)
		{
			Point p0 = _elem->point(0);
			Point p1 = _elem->point(1);
			Point p2 = _elem->point(2);
			Point p01 = (p1 - p0) / ((Real)(_devide));
			Point p02 = (p2 - p0) / ((Real)(_devide));
			Point p12 = (p2 - p1) / ((Real)(_devide));
			Point p0_temp = Point();
			Point p1_temp = Point();
			Point p2_temp = Point();
			Real aera_temp = _elem->volume() / ((Real)(_devide * _devide));
			for (int i = 0; i < _devide; ++i)
			{
				for (int j = 0; j < 2 * i + 1; ++j)
				{
					if (j % 2 == 0)
					{
						p0_temp = p0 + p01 * i + p12 * (j / 2.0);
						p1_temp = p0 + p01 * (i + 1.0) + p12 * (j / 2.0);
						p2_temp = p0 + p01 * (i + 1.0) + p12 * (j / 2.0 + 1.0);
					}
					else
					{
						p0_temp = p0 + p01 * i + p12 * (j / 2.0);
						p1_temp = p0 + p01 * (i + 1.0) + p12 * (j / 2.0 + 1);
						p2_temp = p0 + p01 * i + p12 * (j / 2.0 + 1.0);
					}

					p = (p0_temp + p1_temp + p2_temp) / 3.0;
					_Disintegrate_aera.push_back(aera_temp);
					_Disintegrate_centre.push_back(p);
				}
			}
		}
		else if (_elem->type() == QUAD4)
		{
			Point p0 = _elem->point(0);
			Point p1 = _elem->point(1);
			Point p2 = _elem->point(2);
			Point p3 = _elem->point(3);
			Point p01 = (p1 - p0) / ((Real)(_devide));
			Point p03 = (p3 - p0) / ((Real)(_devide));
			Point p32 = (p2 - p3) / ((Real)(_devide));
			Point p0_temp = Point();
			Point p1_temp = Point();
			Point p2_temp = Point();
			Point p3_temp = Point();
			for (int i = 0; i < _devide; ++i)
			{
				Point p_i = (p3 - p0 + (p32 - p01) * i) / ((Real)(_devide));
				Point p_iplus1 = (p3 - p0 + (p32 - p01) * (i + 1.0)) / ((Real)(_devide));
				for (int j = 0; j < _devide; ++j)
				{
					p0_temp = p0 + p01 * i + p_i * j;
					p1_temp = p0 + p01 * (i + 1.0) + p_iplus1 * j;
					p2_temp = p0 + p01 * (i + 1.0) + p_iplus1 * (j + 1.0);
					p3_temp = p0 + p01 * i + p_i * (j + 1.0);

					p = (p0_temp + p1_temp + p2_temp + p3_temp) / 4.0;
					Real aera_temp = 0.5 * (((p1_temp - p0_temp).cross(p2_temp - p0_temp)).norm() + ((p2_temp - p0_temp).cross(p3_temp - p0_temp)).norm()); 
					_Disintegrate_aera.push_back(aera_temp);
					_Disintegrate_centre.push_back(p);
				}
			}
		}
		else
		{
			mooseError("Too much point，can't disintegrate!");
		}
	}

	else if (dim == 1) //2D
	{
		Real aera_temp = _elem->volume() / ((Real)_devide);
		for (int i = 0; i < _devide; ++i)
		{
			p = _elem->point(0) + (_elem->point(1) - _elem->point(0)) / (2 * _devide + 0.0) * (2 * i + 1.0);
			_Disintegrate_aera.push_back(aera_temp);
			_Disintegrate_centre.push_back(p);
		}
	}

	else // 1D
	{
		//NO NEED
	}
}
